from os.path import join, abspath, basename, splitext
from glob import glob
from datetime import datetime, date, timedelta
from zipfile import ZipFile
import warningsECOSTRESS Workshop Fall 2022
Importing Libraries
These are some built-in Python functions we need for this notebook, including functions for handling filenames and dates.
We’re using the rioxarray package for loading raster data from GeoTIFF files, and we’re importing it as rxr. We’re using the numpy library to handle arrays, and we’re importing it as np. We’re using the rasterstats package for zonal statistics.
!mamba install rasterstats -q -y Package Version Build Channel Size
────────────────────────────────────────────────────────────────────────────────────
Install:
────────────────────────────────────────────────────────────────────────────────────
+ rasterstats 0.17.0 pyhd8ed1ab_0 conda-forge/noarch 18kB
+ simplejson 3.17.6 py39hb9d737c_2 conda-forge/linux-64 107kB
Change:
────────────────────────────────────────────────────────────────────────────────────
- gdal 3.5.0 py39hc691d54_4 installed
+ gdal 3.5.0 py39h85832e7_4 conda-forge/linux-64 2MB
- hdf5 1.12.1 nompi_h2386368_104 installed
+ hdf5 1.12.1 nompi_h4df4325_104 conda-forge/linux-64 4MB
- krb5 1.19.3 h3790be6_0 installed
+ krb5 1.19.3 h08a2579_0 conda-forge/linux-64 2MB
- libgdal 3.5.0 hc0ebe42_4 installed
+ libgdal 3.5.0 h59d0e54_4 conda-forge/linux-64 14MB
- libnghttp2 1.47.0 h727a467_0 installed
+ libnghttp2 1.47.0 hff17c54_1 conda-forge/linux-64 844kB
- libssh2 1.10.0 ha56f1ee_2 installed
+ libssh2 1.10.0 hf14f497_3 conda-forge/linux-64 239kB
- pycurl 7.45.1 py39hd73adbb_2 installed
+ pycurl 7.45.1 py39h9297c8b_2 conda-forge/linux-64 75kB
- python 3.9.13 h9a8a25e_0_cpython installed
+ python 3.9.13 h2660328_0_cpython conda-forge/linux-64 28MB
- tiledb 2.9.5 h1e4a385_0 installed
+ tiledb 2.9.5 h3f4058f_0 conda-forge/linux-64 4MB
Upgrade:
────────────────────────────────────────────────────────────────────────────────────
- ca-certificates 2022.6.15 ha878542_0 installed
+ ca-certificates 2022.9.24 ha878542_0 conda-forge/linux-64 154kB
- certifi 2022.6.15 py39hf3d152e_0 installed
+ certifi 2022.9.24 pyhd8ed1ab_0 conda-forge/noarch 159kB
- cryptography 37.0.2 py39hd97740a_0 installed
+ cryptography 38.0.3 py39h3ccb8fc_0 conda-forge/linux-64 2MB
- curl 7.83.1 h7bff187_0 installed
+ curl 7.85.0 h2283fc2_0 conda-forge/linux-64 93kB
- libcurl 7.83.1 h7bff187_0 installed
+ libcurl 7.85.0 h2283fc2_0 conda-forge/linux-64 359kB
- libpq 14.4 hd77ab85_0 installed
+ libpq 14.5 he2d8382_0 conda-forge/linux-64 3MB
- libzip 1.8.0 h4de3113_1 installed
+ libzip 1.9.2 hc929e4a_1 conda-forge/linux-64 99kB
- openssl 1.1.1p h166bdaf_0 installed
+ openssl 3.0.7 h166bdaf_0 conda-forge/linux-64 3MB
- postgresql 14.4 hfdbbde3_0 installed
+ postgresql 14.5 ha7cec9f_0 conda-forge/linux-64 6MB
Summary:
Install: 2 packages
Change: 9 packages
Upgrade: 9 packages
Total download: 69MB
────────────────────────────────────────────────────────────────────────────────────
Preparing transaction: ...working... done
Verifying transaction: ...working... done
Executing transaction: ...working... doneimport rioxarray as rxr
import numpy as np
from rasterstats import zonal_statsWe’re using the geopandas library to load vector data from GeoJSON files, and we’re importing it as gpd. We’re using the shapely library to handle vector data and the pyproj library to handle projections.
import geopandas as gpd
from shapely.geometry import Point, box
from shapely.ops import transform
from pyproj import TransformerWe’re using the pandas library to handle tables, and we’re importing it as pd.
import pandas as pdWe’re using the seaborn library to produce our graphs, and we’re importing it as sns. We’re using the hvplot library to produce our maps. We’re using the matplotlib library to handle plotting figures, and we’re importing it as plt.
import seaborn as sns
import hvplot.xarray
import hvplot.pandas
import matplotlib.pyplot as pltDefining Constants
These constants define the dimensions of our figures. Feel free to adjust these to fit your display.
FIG_WIDTH_PX = 1080
FIG_HEIGHT_PX = 720
FIG_WIDTH_IN = 16
FIG_HEIGHT_IN = 9
FIG_ALPHA = 0.7
BASEMAP = "ESRI"
SEABORN_STYLE = "whitegrid"
sns.set_style(SEABORN_STYLE)This is the location of the example ECOSTRESS Collection 2 product files.
DATA_DIRECTORY = "/home/jovyan/shared/2022-fall-ecostress-workshop"
print(f"data directory: {DATA_DIRECTORY}")data directory: /home/jovyan/shared/2022-fall-ecostress-workshopLoading an ECOSTRESS Collection 2 Granule
First, let’s trying opening a data layer from a product file.
filename = DATA_DIRECTORY + "/ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01.zip"
print(f"example L2T LSTE file: {filename}")example L2T LSTE file: /home/jovyan/shared/2022-fall-ecostress-workshop/ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01.zipThe granule ID for this granule can be parsed from the filename by dropping the .zip extension.
granule_ID = splitext(basename(filename))[0]
print(f"granule ID: {granule_ID}")granule ID: ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01This product bundle, stored in zip format, contains a number of files, including raster data layers in GeoTIFF format as .tif files, and GeoJPEG browse images as .jpeg files. The GeoTIFF files can be loaded into GIS software, such as QGIS and ArcGIS. The GeoJPEG files can be loaded into Google Earth.
with ZipFile(filename) as zip_file:
for internal_file in zip_file.filelist:
print(internal_file.filename)BadZipFile: File is not a zip fileThe ECOSTRESS Collection 2 tiled products include metadata in JSON format as a .json text file.
with ZipFile(filename) as zip_file:
metadata = zip_file.read(f"{granule_ID}/{granule_ID}.json").decode()
print(metadata){
"StandardMetadata": {
"AncillaryInputPointer": "AncillaryNWP",
"AutomaticQualityFlag": "PASS",
"AutomaticQualityFlagExplanation": "Image matching performed to correct orbit ephemeris/attitude",
"BuildID": "0700",
"CRS": "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +type=crs",
"CampaignShortName": "Primary",
"CollectionLabel": "ECOv002",
"DataFormatType": "COG",
"DayNightFlag": "Day",
"EastBoundingCoordinate": -119.06581640006657,
"FieldOfViewObstruction": "unknown",
"ImageLineSpacing": 70.0,
"ImageLines": 1568,
"ImagePixelSpacing": 70.0,
"ImagePixels": 1568,
"InputPointer": "ECOv002_L2_LSTE_22387_011_20220617T215819_0700_01.h5,ECOv002_L2_CLOUD_22387_011_20220617T215819_0700_01.h5,ECOSTRESS_L1B_GEO_22387_011_20220617T215819_0601_01.h5,ECOSTRESS_L1B_RAD_22387_011_20220617T215819_0601_01.h5",
"InstrumentShortName": "ECOSTRESS",
"LocalGranuleID": "ECOv002_L3T_JET_22387_011_11SKU_20220617T215819_0700_01.zip",
"LongName": "ECOSTRESS Tiled Evapotranspiration Ensemble Instantaneous and Daytime L3 Global 70 m",
"NorthBoundingCoordinate": 35.22502010920582,
"PGEName": "L3T_L4T_JET",
"PGEVersion": "v1.4.8",
"PlatformLongName": "ISS",
"PlatformShortName": "ISS",
"PlatformType": "Spacecraft",
"ProcessingEnvironment": "Linux ceb79f540da6 3.10.0-1160.66.1.el7.x86_64 #1 SMP Wed May 18 16:02:34 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux",
"ProcessingLevelDescription": "Level 3 Tiled Evapotranspiration Ensemble",
"ProcessingLevelID": "L3T",
"ProducerAgency": "JPL",
"ProducerInstitution": "Caltech",
"ProductionDateTime": "2022-10-25T21:39:28.975Z",
"ProductionLocation": "ECOSTRESS Science Computing Facility",
"RangeBeginningDate": "2022-06-17",
"RangeBeginningTime": "21:58:19.929775",
"RangeEndingDate": "2022-06-17",
"RangeEndingTime": "21:59:11.900024",
"RegionID": "11SKU",
"SISName": "Level 3/4 JET Product Specification Document",
"SISVersion": "Preliminary",
"SceneBoundaryLatLonWKT": "POLYGON ((-122.91679849347318 36.98583577347708, -120.16328012620973 40.20637371131396, -117.02446217331683 38.04420848672155, -119.81948100354185 34.92253448622397, -122.91679849347318 36.98583577347708))",
"SceneID": "011",
"ShortName": "ECO_L3T_JET",
"SouthBoundingCoordinate": 34.210027193369996,
"StartOrbitNumber": "22387",
"StopOrbitNumber": "22387",
"WestBoundingCoordinate": -120.29518354417259
},
"ProductMetadata": {
"BandSpecification": [
0.0,
0.0,
8.779999732971191,
0.0,
10.520000457763672,
12.0
],
"NumberOfBands": 3,
"OrbitCorrectionPerformed": "True",
"QAPercentCloudCover": 2.0684334001457727,
"QAPercentGoodQuality": 2.2481680810079174,
"AncillaryNWP": "GEOS.fp.asm.inst3_2d_asm_Nx.20220617_2100.V01.nc4,GEOS.fp.asm.inst3_2d_asm_Nx.20220618_0000.V01.nc4,GEOS.fp.asm.tavg1_2d_lnd_Nx.20220617_2130.V01.nc4,GEOS.fp.asm.tavg1_2d_lnd_Nx.20220617_2230.V01.nc4,GEOS.fp.asm.tavg1_2d_rad_Nx.20220617_2130.V01.nc4,GEOS.fp.asm.tavg1_2d_rad_Nx.20220617_2230.V01.nc4,GEOS.fp.asm.tavg1_2d_slv_Nx.20220617_2130.V01.nc4,GEOS.fp.asm.tavg1_2d_slv_Nx.20220617_2230.V01.nc4,GEOS.fp.asm.tavg3_2d_aer_Nx.20220617_1930.V01.nc4,GEOS.fp.asm.tavg3_2d_aer_Nx.20220617_2230.V01.nc4,GEOS.fp.asm.tavg3_2d_chm_Nx.20220617_1930.V01.nc4,GEOS.fp.asm.tavg3_2d_chm_Nx.20220617_2230.V01.nc4"
}
}To open the temperature layer of this file, we’ll form a Universion Resource Identifier (URI) with the pattern: zip://{filename}!/{granule_ID}/{granule_ID}_{variable}.tif
URI = f"zip://{abspath(filename)}!/{granule_ID}/{granule_ID}_LST.tif"
print(f"URI: {URI}")URI: zip:///Users/halverso/Desktop/2022-Fall-ECOSTRESS-Cloud-Workshop/data/Carpinteria ECOSTRESS Collection 2/ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01.zip!/ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01/ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01_LST.tifWe’re using rioxarray to open the surface temperature layer from the L2T_LSTE product on the 11SKU tile covering the Carpinteria Salt Marsh. We’re passing the URI pointing to the GeoTIFF file contained within this zip file. If you unzip this product bundle or download the GeoTIFF file on its own, you can pass the filename of the GeoTIFF file directly into rioxarray.
ST_K_raster = rxr.open_rasterio(URI).squeeze('band', drop=True)
ST_K_raster<xarray.DataArray (y: 1568, x: 1568)>
[2458624 values with dtype=float32]
Coordinates:
* x (x) float64 2e+05 2.001e+05 2.002e+05 ... 3.096e+05 3.097e+05
* y (y) float64 3.9e+06 3.9e+06 3.9e+06 ... 3.79e+06 3.79e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
_FillValue: nan
scale_factor: 1.0
add_offset: 0.0This xarray.DataArray object contains both an array of image values and spatial metadata. The rioxarray package extends xarray with a .rio attribute containing the metadata. Here we’re examining the coordinate reference system (CRS) of this image in the rioxarray metadata.
CRS = ST_K_raster.rio.crs
CRS.to_wkt()'PROJCS["WGS 84 / UTM zone 11N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32611"]]'This image is projected in UTM zone 11 north. Distances in this projection system are in meters. All tiled ECOSTRESS products are projected in UTM, but tiles are projected into different UTM zones, depending on where they are in the world.
We can also check the spatial resolution of the grid cells in this image in the .rio metadata.
cell_width, cell_height = ST_K_raster.rio.resolution()
print(f"cell width: {cell_width} cell height: {cell_height}")cell width: 70.0 cell height: -70.0This image is projected with 70m square pixels, as are all tiled ECOSTRESS products.
To know the observation time for this granule, we’re parsing the timestamp from the filename. This timestamp is given as Coordinated Universal Time (UTC).
datetime_UTC = datetime.strptime(basename(filename).split("_")[-3], "%Y%m%dT%H%M%S")
print(f"date/time UTC: {datetime_UTC:%Y-%m-%d %H:%M:%S}")date/time UTC: 2020-08-26 19:14:53We want to know the centroid coordinate of this tile so that we can adjust the UTC time given to solar apparent time. We’re calculating the centroid as the average of x coordinate values and average of y coordinate values.
centroid_UTM = Point(np.nanmean(ST_K_raster.x), np.nanmean(ST_K_raster.y))
print(f"centroid UTM: {centroid_UTM}")centroid UTM: POINT (254860 3845120)This centroid coordinate is in meters. We want to convert these projected x and y values in meters to geographic latitude and longitude in degress.
centroid_latlon = transform(Transformer.from_crs(CRS, "EPSG:4326", always_xy=True).transform, centroid_UTM)
print(f"centroid lat/lon: {centroid_latlon.wkt}")centroid lat/lon: POINT (-119.67694453530225 34.71877480226906)We’re shifting the UTC time to local time according to longitude.
datetime_solar = datetime_UTC + timedelta(hours=(np.radians(centroid_latlon.x) / np.pi * 12))
print(f"date/time solar: {datetime_solar:%Y-%m-%d %H:%M:%S}")date/time solar: 2020-08-26 11:16:10The hvplot package extends xarray to allow us to plot maps. We’re reprojecting the raster geographic projection EPSG 4326 to overlay on the basemap with a latitude and longitude graticule. We’re using the jet color scheme to render temperature with a rainbow of colors with red meaning hot and blue meaning cool. We’re setting the alpha to make the raster semi-transparent on top of the basemap. We’re filtering out values lower than the 2% percentile and higher than the 98% percentile to make the variation in the image more visible.
ST_CMAP = "jet"
ST_K_map = ST_K_raster.rio.reproject("EPSG:4326").hvplot.image(
geo=True,
cmap=ST_CMAP,
tiles=BASEMAP,
alpha=FIG_ALPHA,
width=FIG_WIDTH_PX,
height=FIG_HEIGHT_PX,
clim=(ST_K_raster.quantile(0.02), ST_K_raster.quantile(0.98)),
title=f"{granule_ID} Surface Temperature (Kelvin)"
)
ST_K_map = ST_K_map.options(xlabel="Longitude", ylabel="Latitude")
ST_K_mapUnable to display output for mime type(s):